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FOREWORD 

The Study of the KaminUtiqula River, originally planned as a 
waste assimilation capacity investigation in 1985, was subsequently 
expanded and included as a Municipal Industrial Strategy for Abatement 
(MISA) pilot site. Inclusion of the MISA objectives for the site 
study expanded the range of Investigation from traditional nutrient 
and oxygen consuming waste concerns to include all known and suspected 
contaminants from point source dischargers to the river. 

Part 1 of the Kaminlstiquia River water quality study presented 
the findings of water quality surveys carried out in 1986 as they 
relate to the assimilative capacity of the lower river. The findings 
focussed on the Impact of oxygen consuming wastes. A traditional 
one-dimensional riverine model was utilized to define the dissolved 
oxygen depletion. 

This report presents the findings of a joint study between the 
Ontario Ministry of the Environment and Environment Canada. The 
thermal structure and hydrodynamics of the lower river were evaluated 
utilizing a dynamic reservoir simulation model. Temperature profiles 
and flow characteristics were predicted as part of the development of 
water quality models for the MISA program. 

An investigation of the impact of toxic organic and inorganic 
wastes 1s underway utilizing the estuary models and will be presented 
in subsequent technical reports. 

In addition to the MISA study activities, the entire Thunder Bay 
near shore area is under investigation as part of the Remedial Action 
Plan (RAP) process. 
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ABSTRACT 

The one-dimens1onal dynamic reservoir simulation model, DYRESM, 
was modified and applied to the Kaministiqula River for predicting 
temperature profiles and flow characteristics as part of the 
development of water quality models for the MISA program. DYRESM, 
which normally simulates dally vertical lake-wide average temperature 
profiles accounting for meteorological conditions, inflow/outflow and 
diffusion, was modified to simulate a number of adjoining river 
segments. The river experiences extensive stratification due to the 
intrusion of cooler Lake Superior water thus necessitating the use of 
a thermal stratification model. Major modifications are on the 
advection and diffusion of heat along the river reaches and the 
capability to simulate thermal stratificiation at branching channels. 

The data for model calibration and verification are obtained from 
an intensive study by the MOE during August 11 to 16, 1986. 
Modification of estimated velocity profiles and wind effects were 
required during model calibration. A five-day confirmation test from 
August 11 to 15, 1986 resulted in good agreement between simulated and 
observed temperature profiles on August 15 and daily surface water 
temperatures at three locations from August 11 to 14. 
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RESUME 

Le models dynamique unidimensionnel de simulation d'un 
reservoir (DYRESM) a ete modifie et applique a la riviere 
Kaministiquia de maniere a pouvoir prevoir les profils des 
temperatures et les caracteristiques du debit, dans le cadre de la 
raise au point de modeles de qualite de I'eau pour la SMID. 
DYRESM, qui normalement simule les profils des temperatures 
moyennes verticales quotidiennes dans 1' ensemble du lac en tenant 
compte des conditions meteorologiques, du debit entrant et du 
debit sortant ainsi que de la diffusion, a 6te modifie afin de 
pouvoir simuler un certain nombre de trongons contigus de la 
riviere. Celle-ci subit une stratification intense, tout comme un 
lac, en raison de 1' intrusion en amont des eaux plus froides du 
lac Superieur; c'est pourquoi il faut se servir d'un niodele de 
lac. Les modifications principales portent sur I'advection et la 
diffusion de la chaleur dans les divers trongons de la riviere et 
sur la capacite de simuler la stratification thennique aux 
chenaux. 

Les donnees pour I'^talonnage et la verification du modele 
proviennent de l'6tude intensive qu'a menee le ministere de 
I'Environneraent du 11 au 16 aout 1986. II a fallu modifier les 
profils de vitesse et les effets du vent au cours de 1 '^talonnage. 
Un essai de verification d'une duree de cinq jours, soit du 11 au 
15 aout 1986, a donne lieu a une bonne compatibilite entre les 
profils simules et les profils observes des temperatures le 15 
aout et les temperatures quotidiennes des eaux de surface a trois 
endroits differents du 11 au 14 aout. 



SIMULATION OF THE THERMAL STRUCTURE OF THE 
KAMINISTIQUIA RIVER, THUNDER BAY, ONTARIO 

1. INTRODUCTION 

The Kam1n1st1quia River located near Thunder Bay, Ontario is 
subject to Industrial pollutant loadings which often cause the river 
water quality, especially dissolved oxygen concentrations, to fall 
below desired levels. Application of a riverine water quality model 
would normally be sufficient to determine viable solutions. However, 
the Kaministlquia River is unusual since cooler and cleaner Lake 
Superior water Intrudes upstream along the river bottom, which creates 
a vertical thermal structure with a distinct thermocline similar to 
that observed in lakes. This phenomenon also results in both a 
horizontal and a vertical gradient of dissolved oxygen concentration 
since the polluted water is warmer and flows downstream nearer the 
surface. A river water quality model that accounts for vertically 
varying concentration profiles is therefore required. 

Prior to simulating water quality constituents, the flow 
characteristics must be known. In this study, the one-dimensional 
dynamic reservoir simulation model, DYRESM, which simulates lake-wide 
average temperature profiles, 1s modified and applied to the 
Kamlnistiqula River to determine the flow characteristics and water 
temperature. The river is divided into a number of connected segments 
which are simulated in turn using DYRESM in six-hour time steps. 
Permitting inflow and outflow at all depths at segment end points 1s 
one of the two major modifications required in the model. The other 
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modification Involves a new model design to simulate branching, 
stratified flows, as the Kaministiquia River contains three branches, 
all flowing Into Lake Superior. 

2. GOALS AND OBJECTIVES 

The overall goal of the Kaministiquia River Water Quality 
Modelling Study Is to develop and verify a water quality model with 
predictive capability for the assessment of possible management 
strategies for the MISA program on pollution control for rivers of the 
same type. The objective Is to develop a DO-BOD model incorporating 
not only the multi-source, heated effluent conditions but also the 
modulations on the DO-BOD concentrations by the intrusion of the 
relatively cooler and denser lake water. 

At this stage of development, we have just completed the testing 
on the physical flow and temperature modules. This report is an 
interim document limited to showing only the results of the physical 
modules. The results are crucial in the subsequent design of the 
DO-BOD model (Lam et al. , 1987). The DO-BOD model results, however, 
will be presented in another report as a sequel to this one. The 
present report describes the major modifications of the DYRESM 
subroutines as well as the results and findings of the water 
temperature simulation. 
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3. DYNAMIC RESERVOIR SIMULATION MODEL (DYRESM) 

The dynamic reservoir simulation model, DYRESM, was developed to 
predict the vertical variation of temperature and salinity in medium 
sized reservoirs (Patterson et al. , 1978), An ice/snow model was 
added to DYRESM by Patterson and Hamblin (1984) for ice and snow 
simulations but is not needed for this study. The inflow and outflow 
subroutines were replaced by a new simple routine which uses inputted 
flow rate profiles to simulate Inflow and outflow for each river 
segment. In this section, a description of DYRESM, which Is taken 
mostly from McCrimmon (1988), and the modifications for application to 
the Kamlnistiquia River are presented. 

The one-dimensionality of DYRESM Is based on a Lagrangian layer 
structure composed of a number of slab-like layers. Each layer is of 
variable thickness and is assumed to be homogeneous throughout. To 
simulate changes In a reservoir the layers can change their volume, 
temperature and vertical location but not their relative position. 
For example, to model inflow, some layers are thickened or Inserted 
and the layers overlying them are shifted upwards because of the 
underlying Increase In volume. Similarly, outflow Is simulated by 
layer depletion or removal and the downwards shifting of overlying 
layers. Layers can also be amalgamated to simulate mixing. The 
original mainline program, DYRESM, Is converted into a new mainline 
program, KAMRIV, and subroutine, SDYRESM. KAMRIV calculates fixed 
parameter values, controls the input, output and six-hour loop, and 
calls subroutine SDYRESM for each river section In turn within the 
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s1x-hour loop. SDYRESM controls the sub-six-hour loop within which 
Inflow and outflow are simulated and routines 1) HEATR, 2) MIXER and 
3) FCT and DIFUSE are called to simulate meteorological forcing, 
surface mixing and diffusion, respectively. Also used are service 
routines DENSITY, RESINT and THICK to calculate the density, surface 
area and volume (or thickness) of layers and to ensure that layer 
volumes stay within prescribed limits. A flow chart of DYRESM Is 
presented in Figure 1. 

3.1 Main Program KAHRIV 

The initial step of KAMRIV is to read in the fixed data, which 
Includes storage, surface area, length and upstream and downstream 
cross-sectional areas and velocities at different depths, the initial 
temperature profiles, simulation and print control parameters, and 
calibrated parameter values. The daily loop is then Initiated by 
reading In dally meteorological data. Following this the six-hourly 
loop begins and then the sub-six-hourly, forcing-mixing-diffusion loop 
is entered for each segment by calls to SDYRESM. Meteorological 
forcing, surface mixing, deep mixing, inflow and outflow are simulated 
with a sub-dally time step between one quarter hour and six hours. 
The time step, which is calculated 1n subroutine HEATR based on the 
rate of transfer of thermal energy at the surface, ls required to 
better simulate mixing time scales. The sub-six-hourly loop begins 
with the meteorological forcing computations performed by HEATR, 
follows with Inflow/outflow, then epilimnetic mixing using MIXER and 
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ends with the turbulent diffusion simulation using FCT and DIFUSE. 
This loop Is repeated for each river section. Then, the next six-hour 
period Is simulated and so on. After the four six-hour periods are 
finished the dally loop Is re-entered. The last function of KAMRIV is 
to print out the results of the day's simulation. A more detailed 
description of the model will be presented below. 

KAMRIV also calculates Interpolation coefficients (used In 
subroutine RESINT) from the Input storage-depth and surface area-depth 
relationships. The service subroutine RESINT is used to calculate 
layer volumes and surface areas from volumes. Therefore, when layer 
volumes or heights are changed by an operation, then RESINT will be 
called. 

After the Initial temperature profile is read by KAMRIV, the 
subroutines DENSITY and THICK are called. DENSITY calculates the 
density of the water In each layer based upon Its temperature using 
the formula 

DENSITY = 0.9998395 + 6.7914*10-5*T - 9.0894*10-6*t2 + 1.0171*10-7*t3 
- 1.2846*10-9*t4 + 1. 1592*10-11*t5 - 5.0125*10-14*t6 (i) 

where T 1s the water temperature {'C) and the density is in units of 
kg/L. Subroutine THICK is used to ensure that layers do not exceed a 
volume, Vniax« ^" order to maintain a reasonable resolution, and are 
not smaller than a volume, V^in* to limit the number of layers 
required. Imberger and Patterson (1981) found that reasonable values 
are Vmin - S/N and V^ax ' 2 Vn,in where S 1s the lake capacity 
and N Is the maximum number of layers of 150. When a layer Is larger 
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than Vniax* ^t is divided into an appropriate number of layers, all 

less than Vmax* When a layer is less than V(nin» ^t is amalgamated 

with its smallest neighbour to form one layer, then rechecked against 
Vmin *n<^ Vniax- 

3.2 Subroutine HEATR 

KAHRIV enters the daily loop after the above calculations and 
reads the dally data. Immediately after this, the six-hourly loop 
begins, subroutine SDYRESM is called, and the sub-six-hourly loop is 
initiated with the calling of subroutine HEATR. HEATR performs the 
meteorological forcing computations, which consist of calculating: 
transfer of heat at the lake surface due to the shortwave radiation 
incident on the surface, evaporative heat flux, conductive heat flux 
and longwave absorption heat flux. The dally data required for this 
subroutine are: 

SW = Incident short wave radiation (KJ/m2/d) 

SRAT = Sunshine ratio (sunshine hours/daylight hours) 

14 « Average daily a1r temperature ('C) 

SVPD = Average dally air vapour pressure (mb) 

U6 - Average daily wind speed (m/s) 

RAIN = Daily total precipitation (mm) 

All of the calculated sources and sinks of heat at the surface 
are proportioned over each time period and then applied over the 
entire day. The daily Input of shortwave radiation, however, Is 
assumed to apply only over the last 12 hours of a day and, therefore. 
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1s divided by 48 and applied equally over each of the 48 quarter hour 
time periods of a half day. The Input of total dally precipitation Is 
divided equally and applied over each of the 96 quarter hour time 
periods In a day. 

The evaporative and conductive heat fluxes are based on the 
fluxes of moisture, E, and heat, H, given as: 

E = - C/L*U6(f - fj)*p (2) 

H - - Ch*U6{T - T^)*p*Cp (3) 

1n which Ch and C^ are the bulk transfer coefficients, L Is the 
latent heat of evaporation, T is temperature, f is the humidity, p Is 
density, Cp 1s the specific heat of water and subscript s refers to 
the surface. These fluxes are based upon bulk aerodynamic formulae. 
The constants Ch and C^ range from 1.3*10-3 to 1.5*10-3. ^ value 
of 1.4*10-3 was assumed by Imberger and Patterson (1981). The 
formulations of equations 2 and 3 are coded in DYRESM as: 

XEV « - 3.9*U6*(SVP0 - SVPD) (W/m2) (4) 

XCO -= R*XEV = 2.535*U6*(Tg - T4) (W/m2) (5) 

where XEV Is the evaporative heat flux, XCO is the conductive heat 
flux, SVPO Is the average daily saturation vapour pressure at Tg and 
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R 1s the Bowen ratio. DYRESM assumes a positive heat flux 1s directed 
towards the lake surface. 

Longwave radiation heat fluxes are both emitted from the water 
surface and absorbed at the water surface. The emitted radiation Is: 

QW » - 5.53*10-8*(T^ + 273)4 ^gj 

where QW Is the emitted longwave radiation (W/m^). The longwave 

radiation absorbed In the water, QAC (W/m2), is the sum of the 

radiation emitted by clouds plus a fraction of the emitted radiation 

being re-radiated back to the earth and Is modelled as: 

QAC «= (1 + .17C2) (273.1 + 14)6*5.18*10-13 (7) 

where C Is the fraction of cloud cover. 

The penetration of shortwave radiation within the water column Is 
modelled after Beers Law 



q(z) - (1-p) SW e-"l^ 



(8) 



where (1-p) is the fraction of SW penetrating below the surface, ni 1s 
the bulk extinction coefficient and z is the depth. 

The time step used In HEATR Inflow/outflow calculations and MIXER 
Is calculated in HEATR. The minimum time step Is set at 900 seconds 
or one quarter hour. The maximum time step is six hours due to the 
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six hour wind speed components used. In order to prevent the 
turbulent velocity scale, w*. and the mixed layer mean velocity, U, 
from becoming too large In subroutine MIXER, the change In surface 
tenperature before MIXER Is called Is limited to 3'C. The conductive, 
evaporative and longwave heat fluxes are each calculated for one 
quarter hour and summed to determine a change 1n heat In the surface 
layer. Then, If appropriate, the heat Increase due to shortwave 
radiation applied over a quarter hour 1s added. Since the input data 
1s 1n the form of a dally average, it 1s assumed that a temperature 
change due to the above heat sources will be the same for each 
subsequent quarter hour period. Therefore, if AT Is the calculated 
temperature change In the surface layer over one quarter hour, then 
the time step Is limited to 



nq(T) = 3.0/aT (9) 



where nq Is the integer number of quarter hour periods with a 
minimum of one. The time step, nq, must also satisfy the mean 
velocity criterion of 



nq(U) - 0.1 h/u*2 (10) 



where h is the mixed layer depth from the previous time step, u* is 
the current wind shear velocity and U is the mixed layer mean 
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velocity. This criterion is used to ensure the mean velocity increase 
Is limited to 0.1 m/s over the previous value. The temperature 
increment for each layer over nq quarter hours Is then calculated 
and added to the existing temperature. 

The surface layer is adjusted at the end of HEATR to account for 
precipitation and evaporation. The precipitation is a dally input and 
the water level change due to evaporation, WLOST, is also calculated 
dally as 

WLOST = XEV/L (11) 

where L is the latent heat evaporation (=2.453*109 J/m^). 
3.3 Inflow/Outflow 

A new, simple Inflow/outflow routine was developed to replace the 
original DYRESM subroutines, which are inappropriate for the 
tcamlnlstiquia River. Inflow/outflow were originally simulated In the 
dally loop, after the sub-daily loop, but, due to the small retention 
time of the river sections (12 hours or less for period of study), the 
new inflow/outflow routine is used after HEATR within the sub-six- 
hourly loop. 

The new inflow/outflow routine is based on one metre thick 
layers, which are created using RESINT prior to the routine. Upstream 
and downstream velocities, U, and cross-sectional areas. A, for each 
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metre thick layer of each river segment are model inputs. Positive 
and negative velocities are used to Indicate a downstream (D/S) or 
upstream (U/S) direction, respectively. After a time step. At, the 
new volume, V, of a layer is calculated as follows 

V(t + At) - {%fs^un ' "d/S^/S + Q ) ^<^ + V(t) (12) 

where Q is the rate of flow from an external source, such as an 
Industrial outfall. Similar to equation (12), a new layer 
temperature, T, due to inflow/outflow is calculated as 

T(t+At) . [Uu/s%/sTu/S-UD/s'^D/sVs'^Q^Q>"*T<t)V(t)]/V(t+At) (13) 

where 

^U/S ■ ^ °^ *^^ ^^^ segment layer for Uy/s > 0. 

T{t) of the current segment layer for Ujj/s < 0. 
^D/S ~ ^^^^ °^ ^^^ current segment layer for Up/s > 0. 

T of the D/S segment layer for Uq/s < 0. 
Tq = T of external flow 

A time step of At * 15 minutes is used since it equals the smallest 
possible At (the time step of the sub-six-hourly loop). Equations 12 
and 13 are applied to each layer and then, if nq Is larger than At, 
the above procedure is repeated for a total of nq/At times. After 
each At time step the layers volume and thickness change and. 
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therefore, after each At step the one metre thick layers are recreated 
using RESINT. 

3.4 Subroutine MIXER 

Subroutine MIXER Is called after the Inflow/outflow to simulate 
mixing In the epillmnlon. If surface cooling occurred In previous 
routines, then the density of the surface layer may have become 
greater than that of underlying layers, causing the density structure 
to become unstable. To maintain a stable structure MIXER performs 
layer amalgamation for the surface layer with less dense underlying 
layers. 

If layer amalgamation is performed, then the centre of gravity 
win be lowered, releasing potential energy. This potential energy 
per unit area Is the buoyancy flux, APE. If this energy flux Is 
sufficient, then further mixing Is performed, which Increases the 
mixed layer thickness. MIXER then combines any residual energy with 
the wind power flux to attempt further mixing. Layer amalgamation is 
performed until the available energy Is less than the energy required 
to mix with the next layer. Any residual energy Is then stored for 
use In the next mixing event. 

If, after the mixing phases, the Interface is too thin, It will 
be unstable to shear. Subroutine KH Is called to account for the 
formation of the Kelvln-Helmholtz billows at the Interface. If the 
billows are large enough, KH forms at least six layers over the shear 
zone and mixes from the interface outwards to form a linear density 
gradient across the billow thickness. 
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The last function of MIXER 1s to amalgamate layers of equal 
density Irrespective of the volume constraints. The purpose of this 
procedure Is to reduce the computations that will be required 1n the 
diffusion calculations that follow MIXER. 

3.5 Diffusion Subroutines 

Using subroutines FCT and DIFUSE, the final step of the 
sub-six-hourly loop simulates the mixing in the hypolimnion by 
turbulent diffusion. The constant flux model 

dTi 1 d_ ,. ^ . dTK ,,,, 

dt = a;^ di (^p1^ ir^ (14) 

Is the form of the diffusion equation being solved where p Is the 
density, A Is the layer surface area, E Is the eddy diffusivity and 
subscript 1 denotes the layer number. The vertical diffusion 
coefficient is calculated as 

E « ai*H2/(T^*S) (15) 

where a^ Is a mixing efficiency constant, H is the total depth, T^ 
is a mixing time scale and S is a stability factor. The diffusion 
constant aj is actually a function of the basin shape, stratification 
and forcing history but through experience (Imberger and Patterson, 
1981) a constant value of 0.048 produces good results provided the 
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basin shape 1s not too contorted. The stability factor Is calculated 
as 

s ■ f-af (16) 



where 



dz - h^^j - h^ 



h^ " height to centre of layer 1 

p^ ■ density of layer 1 

P - difference in density from top to bottom of reservoir 

H - total lake depth 

FCT is used to calculate equation (15) and the right side of 
equation (14) for each layer. DIFUSE is then used to determine the 
redistribution of heat by ensuring no reversals in the temperature 
gradient occur. This is accomplished by amalgamating two layers when 
the slope of T changes sign across the two layers. 

The mixing time scale, T^, employed in FCT when solving 
equation (15) is calculated as 

\ ' E^tPs-^V (17) 

where Pj is the rate of working of the inflowing water, P^, is the 
rate of working of the wind and E is the potential energy locked in 
stratification. 
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Subroutines FCT and DIFUSE are called for each quarter hour 
period 1n the sub-six-hourly time step. The sub-slx-hourly loop is 
then repeated until the entire day has been covered. The final step 
of the daily loop Is to print the results of the daily calculations. 
The daily loop is repeated for the following days until the end of the 
simulation period. 

4. DATA BASE 

Most Of the meteorological, hydrological and llmnological data 
used for model Input or as observations for model calibration and 
verification were supplied by MOE. In this section, the available 
data, calculated and estimated data and assumptions related to the 
data and model are presented. 

The Kaminlstiquia River Is located in northern Ontario near 
Thunder Bay. The stretch of the river under Investigation extends 
from the river outlet at Lake Superior to approximately 10 kilometres 
upstream, Includes the McKellar and Mission River branches, and Is 
depicted in Figure 2. The points A through P in Figure 2 indicate the 
locations at which parameter measurements were made by MOE. By using 
these points and the added point Z, which Is the location of the 
main pollutant source, as boundaries, 16 river segments were created 
for modelling purposes. 

The majority of time-dependent measurements, such as water 
temperatures and velocities, were made between August 11 and 15 of 
1986. This period was, therefore, chosen as the simulation period. 
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The dally meteorological data required by the model for heat flux 
and mixing calculations are: 

wind speed (daily average) (m/s) 

wind velocity component along longitudinal axis (s1x-hour 

average) (m/s) 

air temperature ('C) 

vapour pressure (mb) 

shortwave radiation (KJ/m^/day) 

precipitation (mm) 
The dally wind speed, air temperature and precipitation were 
supplied. The s1x-hour wind velocity components were calculated using 
hourly wind speeds and directions. Vapour pressure, e, was calculated 
using relative humidity (%) , H, and air temperature, T, in the 
following equation: 

e = .338639H [ (.00738T+.8072)8+. 000019(1. 8T+48)+.001316] (18) 

The shortwave radiation was estimated using the measurements at Big 
Trout Lake. 

Hypsometric data for each river reach or segment are also used as 
model Inputs. Reach lengths and cross-sectional depth-river width 
measurements were used to calculate depth-surface area relationships 
of each river segment by averaging the upstream and downstream cross- 
sectional width and multiplying by the reach length. The widths at 
cross sections Z, B, D, I, K, M. and P. which were not measured. 
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were estimated using upstream and downstream measurements. River 
segment depth-storage relationships were calculated using the 
depth-surface area values. Also, depth-cross-sectlonal area 
relationships were calculated using depth-width values. 

Derivation of hypsometric data for the river branch nodes GHJ and 
KLN required special attention. The drogue data and temperature 
observations indicate the nodes act more like two separate segments 
than one completely mixed segment. Therefore, the nodes were divided 
into river segments GH and GJ at the Mission River branch and segments 
KL and KN at the McKellar River branch. As an example, branch GHJ is 
pictured in Figure 3. The actual proportions of GHJ are not used in 
Figure 3 but the general schematics behind the construction of the two 
segments GH and GJ are displayed. In this case, the cross section G 
is divided into two parts G^ and Gj first, and then the segment GH 
is formed between cross sections Gh and H and the segment GJ between 
Gj and J. The thermocline of segment GH is depicted as being lower 
than GJ's but this would affect GJ only for the case when the flow 
moves from the segment GH to the upstream segment FG (Fig. 2) and then 
moves from FG to GJ, or vice versa, from GJ to FG, and then onto GH. 
In other words, GH and GJ are separated by an invisible wall in the 
model which is not true in reality {Fig. 2). There are more elegant 
ways of treating a river conjunction but the present method is simple 
and convenient to implement for a model such as DYRESM. 

The cross-sectional area of the upstream cross sections, 6 and K, 
are then proportioned to each river section based on the proportion of 
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flow going to each section. These proportions are: 1) 42% to GH and 
58% to GJ, and 2) 31% to KL and 69% to KN. 

Cross-sectional velocity profiles 1n conjunction with cross- 
sectional areas are used by SDYRESM to calculate inflow and outflow of 
each river section. Observed velocities were obtained at seven 
locations at the surface using a velocity meter and at six locations 
at 2 m depth, five locations at 4 tn depth, and four locations at 6 m 
depth using drogues. Since velocity profiles could not be derived for 
each cross section based solely on these observed values, temperature 
contours were constructed from the August 15 observations to help 
trace apparent water movement and help interpolate or extrapolate 
velocities. 

After the first velocity estimates, conservation of mass was 
checked for each cross section by comparing the observed flow rate 
with the summed product of velocity and cross-sectional area of each 
layer at the cross section. Velocities were modified if mass was not 
conserved. The observed flow rates are presented next. 

Flow rates in the river are used to estimate the velocity 
profiles for inflow and outflow calculations. The flow rate along the 
river increases approximately 2 (m^/s) at Z and 1 (m^/s) at B due to 
industrial effluent loading. The nearest upstream location of dally 
measured flow rates for the Kaminlstiqula River is at Kakabeka Falls, 
which is approximately 35 km upstream of A. These values may not be 
appropriate for use as the flow rate at A since a small river flows 
into the Kaminlstiqula River 6 km downstream of the falls. However, 
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on August 14, a flow rate was measured at Stanley, which 1s 1 km 

downstream of the tributary. By assuming the flow at Stanley 

represents the flow at A, the flow rates for the five-day simulation 

period, August 11 to 15, were calculated as follows: 

0^(1) - Q^p(I-l)*(19.9/17.6) (19) 

where ^l^^^ ' ^^°^ ^^^^ ^^ ^ 

Q|,p(I-l) = flow rate measured at Kakabeka Falls for 

previous day 
I = day 

Note that the ratio of flow at Stanley, August 14 to Kakabeka Falls, 
August 14 Is 19.9/17.6. The previous days flow at Kakabeka Falls is 
used in equation (19) because the travel time from Kakabeka Falls to A 
is approximately one day. The velocity profiles are scaled using the 
variation of flow calculated by equation (19) to account for the 
changing flow rate while keeping the water level constant. 

Daily water surface elevations at each cross section were 
available from staff guage readings. The water level varied less than 
0.15 m along the stretch of river over the entire simulation period. 
Therefore, to simplify the model, the water level was assumed to be 
constant. An approximate average value of 183.5 m was used. It is 
noted, however, that changes In water level at other times could be 
Important and must be Incorporated In the model. 
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Water temperature observations consisted of 1) temperature 
profiles at the right, middle and left of each cross section except 
for M and 2 on August 15, 2) four-hourly surface temperatures at all 
stations except Z from August 11 to 15, and 3) hourly surface 
temperatures at stations A, B, D, 6, I, K and P from August 11 to 14. 
At each cross section the right, middle and left profiles were 
averaged. The temperature profiles for each river segment were then 
estimated as the average of upstream and downstream cross section 
average profiles (e.g. the temperature for segment GH Is the average 
of those at cross section Gh and H, Fig. 3). These river section 
profiles are used as Initial conditions and for comparison to 
simulations. The profiles at the river boundary cross sections of A, 
I, M and P were assumed to represent the temperatures of the water 
flowing downstream Into section AZ and flowing upstream into segments 
HI, LM and OP, respectively. These Inflow temperature profiles were 
scaled every six hours of the simulation period using six hour 
averages of the hourly surface temperatures at sections A, I and P and 
the assumption that the temperatures at M equal P. This scaling of 
the inflow temperatures is made to better reflect their change with 
time. ^ 

5. RESULTS 

5.1 Model Calibration 

The objective of calibrating KAMRIV is to reproduce the observed 
temperature profiles of each river section. The profiles measured on 
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August 15, 1986 were selected for calibration since such detailed 
observations were available only for this day. By assuming, 
hypothetical ly, that the same Input data of August 15 hold for five 
consecutive days, a steady state simulation was performed. 

The first simulation prior to any parameter modifications to the 
original DYRESM model resulted In reasonably good temperature profiles 
upstream of the Mission River branch. Below the branch, most river 
segments were simulated as having the top 4 to 5 metres completely 
mixed with temperatures equal to approximately the observed surface 
value. Conversely, the observed profiles indicated little mixing 
because the temperature, In general, decreased 4 to 6'C in the top 4 
to 5 metres. Therefore, it appeared too much surface mixing was 
simulated. This overestimation of wind mixing was not surprising, 
because DYRESM was originally intended for lakes where surface mixing 
Is subject to wider exposure to wind than rivers. 

A wind reduction factor was therefore required to decrease the 
surface wind mixing. In fact, reducing the winds to zero greatly 
Improved the simulations of the lower segments, but the computed 
temperatures remained too high (up to 5'C higher than the observed) In 
the top 3 to 4 metres. The upper river segments showed negligible 
changes. Above cross section G, the average detention time is three 
times smaller than that of the lower segments. This Indicates the 
upper segments are controlled more by inflow/outflow, which probably 
explains why their temperature profiles were not affected by changing 
wind values. Observations show that below cross-section G the 
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thentiocline rises relatively quickly. Also, the surface flow rates 
-Increase down the river Indicating the downstream flowing water 1s 
gradually being forced more to the surface. To further calibrate the 
model, the velocity profiles of the lower segments were modified to 
reflect the above findings. As before, the modified velocity profiles 
were subject to satisfying the mass conservation constraint. 

Through modification of the velocity profiles and reduction of 
the wind effects, reasonably good agreement (rms=0.97*C and average 
relative error = 5.7%) between the simulated and observed temperature 
profiles of each segment was achieved. These profiles are plotted In 
Figure 4. Two notlcable exceptions are segments AZ and ZB, which were 
both simulated as virtually completely mixed but the 'observed' 
profiles showed strong stratification. In the model, the segment ZB 
Is completely mixed due to the Introduction of 37°C effluent at the 
rate of 2 m^/s at the bottom layer. This less dense water causes 
upward mixing and because ZB Is a small segment the mixing proceeds to 
the water surface. Due to upstream flow along the river bottom, the 
segment AZ also receives some warm bottom water from the segment ZB, 
which causes AZ to also be completely mixed. 

On the other hand, the 'observed' profiles for segments AZ and ZB 
were constructed without these considerations. In fact, the 
'observed' temperature profile for the added cross section Z needed to 
be defined first and was done by simply averaging the observed 
temperatures at cross sections A and B. Then, the 'observed profile' 
for segment AZ was obtained by averaging those at cross sections A and 
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Z, and subsequently that for ZB by averaging those at Z and B. 
Therefore, the temperature profiles measured at station B has direct 
Influence on the estimated or 'observed' profiles in both segments AZ 
and ZB. Thus, we do not know whether the observed profiles for AZ and 
ZB as shown in Fig. 4 represent reality or not. 

5.2 Model Confirmation 

The modified DYRESM was previously calibrated using the August 15 
Input data and observations. To attempt to confirm the model, the 
Kaministlquia River was simulated from August 11 to 15. Since the 
August 15 observed temperature profiles, which were used in the model 
calibration, are also used here, this cannot be considered a true 
model verification simulation. However, the dally varying meteorolo- 
gical data, flow rates (velocity profiles) and upstream and downstream 
temperature profiles are used, which provides for a good confirmation 
of the model. More comprehensive verification of the model is pending 
on the availability of an independent set of data, e.g. the 1987 data. 

The resulting temperature profiles on August 15 of the confinna- 
tlon simulation are presented 1n Figure 4. The simulations compared 
reasonably well (rms = O.gS'C and average relative error = 5.8%) with 
the observed profiles and were slightly better for some profiles and 
slightly worse for others when compared to the calibration results. 
Recall that the rm5=0.97"C and average relative error 1s 5.7% for 
model calibration. 
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Surface temperatures at cross sections D, G and K were measured 
over August 11 to 14 and were compared to the simulated values to 
further confirm the model. KAMRIV outputs the surface temperature of 
each river segment after each six-hour period. For comparison, 
six-hour averages of the hourly and four hourly observations were 
calculated. To estimate the simulated temperatures at the desired 
cross sections, the temperatures of river segments imnedlately 
upstream and downstream of the cross sections were averaged, as 
described earlier. Therefore, the simulated surface temperature at D 
equals the average of CD and DE, that at G equals the average of FG 
and the average of GH and GJ, and that at K equals the average of JK 
and the average of KL and KN. 

The computed versus simulated surface temperatures are plotted in 
Figure 5. Statistics comparing these values were calculated and are 
presented in Table 1. In general, as seen 1n Figure 5 and by the low 
r values, the computed and simulated temperatures did not correlate 
well. However, the small values of the root mean square error, as 
presented in Table 1, indicate good results with respect to the 
magnitude of the temperatures. These results suggest that DYRESM is 
not well suited for simulating small time scale fluctuations but is 
good for simulating mean magnitudes of temperature over a longer time 
scale. For example, if dally average values are used, then the rms 
value decreases to 1.0 and the r value increases to 0.53 for D, G and 
K combined, which Indicates better agreement. 
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TABLE 1. Computed Versus Observed Surface Temperature Statistics, 
August 11-14 



Cross Section Regression Mean Mean Root Mean 

Coefficient Computed Observed Square Error 



D 


0.45 


20.0 


21.0 


1.18 


G 


-0.03 


20.4 


20.4 


1.06 


K 


-0.02 


19.8 


18.6 


1.54 


D.G and K 


0.23 


20.1 


20.0 


1.23 


Combined 










Daily Averages 


0,53 


20.1 


20.0 


1.04 


D, G and < 











5.3 Sensitivity Analysis 

The effects of wind on temperature profiles was investigated 
during model calibration. To evaluate the relative importance of 
surface heat fluxes, a sensitivity analysis involving shortwave 
radiation was performed. The daily shortwave radiation was set to 
zero and the model was rerun for August 11 to 15. Compared to the 
confirmation results of August 15, using zero shortwave radiation 
resulted in surface temperatures decreasing 0.2 to 1.2''C and bottom 
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temperatures decreasing 0.2 to 0.4'C. Also, similar to the wind 
sensitivity analysis, the lower river segments were most affected by 
the reduced radiation. These results Indicate that the surface heat 
fluxes are a necessary part of the model since simulations 
significantly worsened without shortwave radiation. The daily total 
fluxes from the verification run are presented in Table 2. 
Evaporative and conductive fluxes are not included because they are 
zero when wind speed is zero. 

TABLE 2. Dally Surface Heat Fluxes (KJ/m^/day) 



Date 


Shortwave4. 


Longwave t 


Longwave^ 


Net Influx 


Aug. 11 


24250. 


35101. 


24776. 


+13925. 


Aug. 12 


17460. 


35102. 


25214. 


+7572. 


Aug. 13 


12610. 


35007. 


27964. 


+5567. 


Aug. 14 


7760. 


34983. 


29963. 


+2741. 


Aug. 15 


14550. 


34929. 


28208. 


+7829. 



Downward shortwave radiation 1s reduced 3% to account for surface 
reflection and It Is this reduced value that is presented In Table 2. 

A sensitivity analysis for the water extinction coefficient was 
performed. It was found that doubling the coefficient resulted In 
temperature changes less than I'C at the surface. This Indicates a 
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relative insens1t1v1ty to surface heat flux compared to Inflow/outflow 
effects. 

6. CONCLUSION 

DYRESM was modified for application to the Kamlnlstlqula River. 
The main modification Involved the replacement of the inflow and 
outflow routines by a simple routine, which utilizes inputted velocity 
and cross-sectional area profiles to calculate inflow/outflow at the 
upstream and downstream boundaries of each river segment. DYRESM was 
used to simulate all segments, including branching segments. In turn 
every six hours. 

Calibration of the model involved reducing the wind effect to 
zero and modifying velocity profiles to reflect a rising thermocline 
and increasing surface velocities with distance downstream. Good 
agreement between simulated and observed temperature profiles was 
achieved. The wind reduction improved the simulation of segments 
below the first river branch by reducing the surface wind mixing. 
Upstream of this branch the simulations did not change notlcably as 
the flow appears to dominate In these segments, which have three times 
smaller detention times. However, it is possible that the winds 
should have been reduced only a certain fraction {due to shelter 
effects of surrounding terrain) and then the velocity profiles should 
have been further calibrated. Also, the existing mixing formulations 
In DYRESM might be inappropriate and may require modification or 
calibration to allow for wind effects. 
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The river was simulated from August 11 to 15 for model 
confirmation. For comparison purposes, the observed temperature 
profiles of each station for August 15 and the observed hourly and 
four-hourly surface temperatures at cross sections D, G and IC were 
used. Good agreement for the profiles was achieved. On a six-hour 
average basis the surface temperature did not correlate well but the 
root mean square errors were approximately only 5% of the mean 
temperature, which is reasonably good. On a daily average basis the 
surface temperatures agreed better, due probably to the fact that 
DYRESM was originally written for daily simulations. For example, 
heating due to shortwave radiation Is simulated In DYRESM over the 
last 12 hours of a day but In reality shortwave heating occurs during 
the middle of the day. Therefore, reasonable correlations over time 
steps of less than one day cannot be expected. 

The modified DYRESM model was found to be practical for 
estimating the Kamlnlstiquia River flow characteristics and, 
consequently, for simulating water temperatures. The model was 
calibrated and confirmed for the period August 11-15, 1986, which was 
during the lowest flow period of 1986. To apply this model during 
another period, which will likely have higher flow rates and water 
levels, new estimated velocity profiles will likely be required as 
well as modification of the model to permit fluctuating water levels. 

On the whole, we believe that the modified DYRESM model now 
provides not only reliable estimates for water temperature and flow 
for use In developing the DO-BOD model, but also the knowledge of the 
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thennal structures along the river during the simulation period that 
are necessary for the design of the spatial compartments and the 
choice of time steps 1n the water quality model. 
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Figure 4. River Segment Temperature Profiles 
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